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Abstract 

We propose a novel numerical method for solving a quadratic vector 
equation arising in Markovian Binary Trees. The numerical method con- 
sists in a fixed point iteration, expressed by means of the Perron vectors 
of a sequence of nonnegative matrices. A theoretical convergence analysis 
is performed. The proposed method outperforms the existing methods for 
close-to-critical problems. 

1 Introduction 

In this paper we study the quadratic vector equation 

x — a + B{x($x), (1) 

where a £ E", B £ E"^" have nonnegative entries, the symbol €5 denotes 
the Kronecker product, and the unknown x is an n-dimensional vector. The 
coefficients a and B are such that the vector e = (1, 1. . . . , 1)-'" is a solution of 

Equation M arises in the study of Markovian Binary Trees (MBT), which 
are a particular family of branching processes used to model the growth of pop- 
ulations consisting of several types of individuals, who may produce offsprings 
during their lifetime. MBTs have applications in biology, epidemiology and also 
in telecommunication systems. We refer to [US] for definitions, properties and 
applications. 

One important issue related to MBTs is the computation of the extinction 
probability of the population, which is the minimal nonnegative solution x* G 
E" of the quadratic vector equation ([T]). 
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The MBT is called subcritical, supercritical or critical if the spectral radius 
p{R) of the matrix R = B{e® I + 1 ® e) is strictly less than one, strictly greater 
than one, or equal to one, respectively. In the subcritical and critical cases the 
minimal nonnegative solution is the vector of all ones, while in the supercritical 
case X* < e, X* ^ e (see Hj and |T]). Thus, only the supercritical case is of 
interest for the computation of x* . 

Several numerical methods have been proposed for computing the vector x* . 
In [2] the authors propose two linearly convergent algorithms, called depth and 
order algorithms. The thicknesses algorithm, still linearly convergent, is pro- 
posed in [3]. In [3] the authors apply the Newton method, which has quadratic 
convergence. A modification of Newton's method has been proposed in [5]. All 
these methods have a probabilistic interpretation, and each of them provides a 
sequence {xk}k of nonnegative vectors, with = (0, . . . ,0)"^, which converges 
monotonically to the minimal nonnegative solution x* . A common feature of all 
these methods is that their convergence speed slows down when the problem, 
while being supercritical, is close to critical, i.e., the spectral radius of R is close 
to one and x* ~ e. Moreover, the accuracy of the approximation deteriorates. 

In this paper we write equation ([T]) in the form x = a+b{x, x) where 6(m, v) := 
B{uiSi v) is the bilinear form defined by the matrix B. If we set y = e — x, the 
latter equation becomes 

y = b{y,e) + b{e,y) -b{y,y). (2) 

The sought solution y* of ([2| is y* = e — x*, where x* is the minimal nonnegative 
solution of ([!]). In the probability interpretation of Markovian Binary Trees, 
since x* is the extinction probability, then y* ^ e — x* is the survival probability. 

Applying a functional iteration directly to ([2]), like Newton's method, gives 
nothing new, since ([2|) differs from ([T]) by a linear change of variable. However, 
the new equation ([2j) can be rewritten as 

y = Hyy, (3) 

where Hy := 6(-, e) + b{e — y, ■). The matrix Hy is nonnegative and irreducible 
if y < e. In particular the solution y* is such that p{Hy ) = 1 and y* is the 
Perron vector of the matrix Hy* . 

This interpretation allows to design a new algorithm for computing y* . To 
this purpose, define the map PV(A/) as the Perron vector of a nonnegative 
irreducible matrix M, so that we may rewrite ([s]) as 

y^VY{Hy). (4) 

The idea is to apply a fixed-point iteration to solve Q, thus generating a se- 
quence {yk}k of positive vectors such that yu+i = P^{Hyk) ^^^d yk converges to 
y* . A suitable normalization of the Perron vector, consistent with the solution, 
is needed to obtain a well-posed iteration. In this way we obtain a new iter- 
ative scheme, which is completely different from classical functional iterations. 
Indeed, the proposed algorithm, unlike known methods, fully exploits the fact 
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that the solution x — e oi the equation ([2]) is known. Moreover, the fixed-point 
iteration at the basis of our algorithm relies on a new interpretation of the solu- 
tion y* in terms of the Perron vector. These differences with respect to classical 
methods lead to great improvements in the numerical solution of MBTs which 
are close to critical. 

We perform a convergence analysis of the fixed point iteration yk+i = 
PV{Hyi^), by giving an expression to the Jacobian of the map y — )■ PV(iJy). 
This expression allows to derive a local convergence result. Moreover, most 
importantly, we prove that, although the convergence of the method is linear, 
the speed of convergence increases as the problem gets close to critical. In the 
limit case of a critical problem, the convergence becomes superlinear. This nice 
behavior is opposite to the one of Newton's method, whose speed of convergence 
is sublinear in the supercritical case, and becomes linear in the critical case. 

A wide numerical experimentation confirms our theoretical analysis. For far- 
from-critical problems the standard techniques are preferable, while for close- 
to-critical problems our method outperforms the existing ones. 

The paper is organized as follows. In Section [2] we state our assumptions on 
the problem. In Section[3]we rewrite the vector equation in terms of an equation 
for the vector y = e — x and discuss the properties of the equation obtained in 
this way. The new algorithm, based on a Perron iteration, is presented in Sec- 
tion |4] The theoretical convergence analysis is performed in Section [5j Finally, 
in Section [6] we present the results of the numerical experiments. Conclusions 
and open issues are addressed in Section [Tj 

2 Assumptions on the problem 

Let a £ M", B G M"^" have nonnegative entries, and consider the quadratic 
vector equation ([T]) where it is assumed that the vector e = (1, 1, . . . , 1)-^ is a 
solution. Let x* € E" be the minimal nonnegative solution of ([I]), i.e., x* < x 
for any other nonnegative solution, where the semi-ordering is component-wise. 
A unique solution x* exists, according to the results of pj Section V.3]. 
We assume that p{R) > 1, where 

R = B{e®I + I®e). 

Under this assumption x* < e, x* ^ e (see [1] and [T]). It is worth pointing out 
that if p{R) = 1, then x* — e, therefore as p{R) is greater than 1 and gets closer 
to 1, then X* approaches to the vector of all ones. 
We introduce the bilinear operator 

: M" X M" M" 

defined as 

b{u, v) — B{u ® v) 

and rewrite ([T]) as 

X = a + b{x, x). (5) 



3 



We assume that for the minimal solution x* of ^ it holds < x* < e, and that 
the Jacobian of the map x ^ x — a — b{x, x) at x* , i.e., / — b{x* , •) — b{-, x*), is 
a nonsingular irreducible M-matrix. Since irreducibility is only determined by 
the nonnegativity pattern of 5(-, •), the irreducibility condition is equivalent to 
requiring that b{e, ■) + b{-,e) is irreducible. Notice that the latter is just another 
notation to represent the matrix R. 

Moreover, we may assume that e^b{e — x*,e — x*) > 0, otherwise b{e — 
X* ,e — X*) = (since B > 0), and the problem is trivial since it becomes a 
linear problem. 

3 The optimistic equation 

A property of equation ([5| which has not been exploited so far in the existing 
literature is that a; = e is a solution. If we set a; = e — y, by using the bilinearity 
of the operator b{-, •) and the property that e = a + b{e, e), equation ^ can be 
rewritten as 

y ^b{y,e) + b{e,y) -b{y,y). (6) 

The trivial solution is ?/ = 0, which corresponds to x = e. We are interested in 
the nontrivial solution < y* < e, which gives the sought solution x* = e — y* . 

In the probabilistic interpretation of Markovian Binary Trees, x* is the ex- 
tinction probability, thus y* = e — x* is the survival probability, i.e., y* is the 
probability that a colony starting from a single individual in state i does not be- 
come extinct in a finite time. For this reason, we refer to ^ as to the optimistic 
equation. 

Notice that ([6| admits the following probabilistic interpretation. The term 
b{y, e) represents the probability that the original individual Ai (for "mother" ) 
spawns an offspring (for "first-born" ) , and after that the colony generated by 
the further offsprings of A^, excluding J", survives. The term b{e,y) represents 
the probability that M spawns T, and the colony generated by survives. The 
term b{y,y) represents the probability that spawns and after that both 
their colonies survive. Thus (|6| follows by the well-known inclusion- exclusion 
principle 

P[A^ or T survives] — ¥[A4 survives] +P[J^ survives]— P [both Ai and survive], 

where V[X] denotes the probability of the event X. 
Equation ^ can be rewritten as 

y = Hyy (7) 

where 

i/, = 6(-,e)+5(e, •)-%,•)■ (8) 

Notice that Hy is the sum of a fixed matrix and a matrix that depends linearly on 
y. Therefore the quadratic operator on the right-hand side of Q is "factored" 
as the product of a matrix which depends on y, and y. 
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An important property is that Hy is a nonnegative irreducible matrix, when- 
ever y < e. Therefore, by the Perron- Frobenius theorem [7], if y < e, has 
a positive eigenvalue Ay = p{Hy), the so-called Perron value, and to Xy corre- 
sponds a positive eigenvector Wy, unique up to a multiplicative constant, the 
so-called Perron vector, so that HyWy = XyWy. Therefore the sought solution 
y* can be interpreted as the vector < y* < e such that p{Hy* ) — 1 and y* is 
a Perron vector of Hy . 

It is worth pointing out that this interpretation of y* in terms of the Perron 
vector allows to keep away from the trivial solution y ~ of Q, since the 
Perron vector has strictly positive elements. 

The formulation of the quadratic vector equation in terms of the Perron 
vector allows to design a new algorithm for its solution. 



4 The Perron iteration 

If we set up a fixed-point iteration or a Newton method for y based on ([6]), 
we get the traditional fixed-point iterations and Newton methods for MBTs [J , 
since what we have done is simply a linear change of variables. Instead, we 
exploit the fact that y* is a Perron vector of the nonnegative irreducible matrix 
Hy. (compare Q). 

To this purpose we introduce the operator 

u = PV(X) 

which returns the Perron vector u of the irreducible nonnegative matrix X. 

Thus, we can devise a fixed-point iteration to compute the solution y* by 
defining the sequence of vectors 

2/fe+i = PV(i/yJ, fc- 0,1,2,..., (9) 

starting from an initial approximation yo. In order to define uniquely the se- 
quence {yk}, we need to impose a normalization for the Perron vector, which 
is uniquely defined up to a multiplicative constant. A possible choice for the 
normalization is imposing that the residual of (|6| is orthogonal to a suitable 
vector w e M", i.e., 

w^{yk+i - b{yk+i,e) - b{e, yu+i) -f b{yk+i,yk+i)) = 0. (10) 

Clearly, this normalization is consistent with the solution of (|6|. We choose 
w as the left Perron vector of the matrix 6(-, e) -I- h{e, •); the rationale for this 
choice is discussed in Section [H 

Given a Perron vector u of Hy^ , the equation to compute the normalization 
factor a such that yk+i — ceu satisfies ( [To| reduces to 

aw u — aw b{u,e) + aw 6(e, u) — a w b[u,u). 
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whose only non-zero solution is 



w'^ {u — b{u, e) — b{e, u)) 
w'^b{u, u) 

Notice that the solution a — corresponds to the trivial solution y ^ Q [x — e)^ 
which we want to avoid. 

The PV(-) operator is defined on the set of irreducible nonnegative matrices. 
If y < e, then the matrix Hy is nonnegative irreducible, therefore the sequence 
yk generated by (|9| is well defined if yu < e for any k. 

In Section [5] we show that the iteration ^ is locally convergent. Therefore, 
if yo is quite close to y* , one can expect that y^ < e for any k. In the case 
where Hy^, is not a nonnegative irreducible matrix, we can define yk+i as an 
eigenvector corresponding to the eigenvalue of Hy^, having maximal real part. 
We call maximal eigenvector this eigenvector. Clearly if Hy^ is a nonnegative 
irreducible matrix, the maximal eigenvector is the Perron vector. We see in 
Section |6] that this concern is not necessary in practice. 

As a starting approximation yo we may choose the null vector. For close 
to critical problems, where y* is close to zero, this choice should guarantee the 
convergence, according to the results of Section [Sj 

The resulting iterative process is summarized in Algorithm 1. 



Algorithm 1 The Perron iteration 
Set fc ^ 
Set yo ^ 

Set w the Perron vector of b{e, ■) + b{-, e) 
while \\Hy^yk - yk\\i > e do 

Set u ^ the maximal eigenvector of Hy^^ 

Compute the normalization factor a — — ("~K".e)-b(e,M)) 

Set yk+i ^ au 

Set fc ^ fc + 1 
end while 
return x = e — yk 



5 Convergence analysis of the Perron iteration 

In this section, we show that the Perron iteration ([9| is locally convergent, and 
its convergence is linear. Moreover, the convergence speed gets faster as the 
problem gets closer to critical. 

5.1 Derivatives of eigenvectors 

It is well known [8] that the eigenvalues and eigenvectors of a matrix are ana- 
lytical functions of the matrix entries in a neighborhood of a simple eigenpair. 
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The following formula for an analytical expression of their first derivatives is 
from Meyer and Stewart P, Theorem 1]. 

Theorem 1. Let A = A{z), A = A(z), u = u{z) be a matrix, eigenvalue and 
associated eigenvector depending on a parameter z G C. Let us suppose that 
A(zo) simple and A'{zq), A'(zo), u'{zq) each exist. Let w — ^/{z) be another 
vector such that w'{zq) exists and let a{u, w) be a function whose value is a real 
scalar constant for all z. Let and be the partial gradients of a{-, •) seen 
as a function respectively of its first and second vector argument only. 
If a^u 7^ for z ~ zq, then the derivative u' of u at z — zq is given by 

u' = -nA-XI)*A'u-a-u.'^ 
a" u 

Here X"^ denotes the so-called group inverse of a singular matrix X, i.e., the 
inverse of X in the maximal multiplicative subgroup containing X. We refer 
the reader to the abovementioned paper for more details on group inverses. 

In fact, very little is needed on group inverses, and the formula can be 
modified slightly in order to replace it with the Moore-Penrose pseudoinverse 
X^ , which is a more canonical tool in matrix computations. 

Theorem 2. With the same hypotheses as Theorem^ let v{z) be the left eigen- 
vector of A{z) corresponding to the eigenvalue A(z). If u 7^ for z = Zq, 
then the derivative u' of u at z — Zq is given by 



^nA-XlViA' X'I)u-.-u.'^ _ _ _ ^^^^ 



(11) 



with A' 



v" A'u 



Proof. The proof is a minor modification of the original proof [6] of Theorem [T] 
By differentiating the identity An = Xu we get A'u + An' = X'u + Am', i.e., 

(yl-A/)u' = -(A'-A7)M. (12) 

By left-multiplying everything by w^, and noting that v^ A = Xv^ , we get 
the required expression for the eigenvalue derivative A'. Moreover, since u is a 
simple eigenvector at z = zq, the kernel of {A — XI) is span(u). Thus from ( 12 1 
we can determine u' up to a scalar multiple of u: 

u' ^ -{A- XlYiA' - X'I)u + Su. (13) 

We shall now use the normalization condition a{u,w) — k to determine the 
value of S. By differentiating it, we get 

cr^{u,w)u' + cr^{u,w)w' ^0. (14) 



Plugging (13) into (14) yields a linear equation for S. □ 
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5.2 Jacobian of the Perron iteration 

The Perron iteration is a fixed-point iteration for the function F{y) := PV{Hy), 
where the function u — PV{X) returns the Perron vector u of the nonnegative 
irreducible matrix X, normalized such that w'^{u — Huu) — 0, where w is a 
fixed positive vector. We can use Theorem [2] to compute the Jacobian of this 
map F . 

Theorem 3. Let y be such that Hy is nonnegative and irreducible. Let u = 
F{y), and let v such that v^Hy — , where A = p{Hy). Then the Jacobian 
of the map F at y is 

where 

aj" ~ w'^{L — b{e — u,-) ~ b{-, e — u)). 

Proof. We shall compute first the directional derivative oi F a.t y along the 
direction a. To this purpose, let us set y{z) := y + az, for any z S C, and 
A{z) = Hy. We have 

d 
dz 

Moreover, set 



A\z)^-Hy = -b{a,-). 



a{u, w) — iv^ {u — &(e, u) — b{u, e) + b{u, m)), 

where w{z) = w for each z (so that w' — 0). The partial gradient of cr(-, •) with 
respect to the first argument is af — uF {L — b{e — w, •) — &(•, e — u)). Plugging 
everything into (ITl|), we get 



a( u 

L ^ {A - Xjy A' ^ — / u 

a{uj \ u J 

/-^l (^-A/)t f/-^) A'u 

"4-) [l-'^) i-Ka.u)) 



= n-^\[A-\Ly \L-^]b{-,u)a. 

From this expression for the directional derivative, it is immediate to recognize 
that the Jacobian is p^. □ 
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5.3 Local convergence of the iteration 

The fixed-point iteration yt+i = F{yk) is locally convergent in a neighborhood of 
y* if and only if the spectral radius of JFy- is strictly smaller than 1. First notice 
that it makes sense to compute the Jacobian using ([T5| in a neighborhood of the 
solution?/*. In fact, afy* = w'^{y* -b{e-y* ,y*) -b{y* ,e-y*)) = w'^b{y*,y*) 
and the latter quantity is positive as w > and b{y*,y*) > 0, as stated in 
Section [2] Moreover, since A = 1 is a simple eigenvalue, the left and right 
eigenvectors v — v* and u — y* cannot be orthogonal. 



By evaluating (151 at j/ = y*, we get 



* *T \ / * *T 



^ \ w^b{y*,y*)J \ v*^y* J ^ ^ ' 

where we have set al'^ = w'^ {I - b{e - y* , ■) - b{- , e - y*)) and A = {b{e-y*, •) + 
6(-,e)-/). 

Let us try to understand what happens to the spectral radius p{JFy- ) when 
the problem is close to critical. 

Theorem 4. Let bt{-,-), t £ [0,1] be an analytical one-parameter family of 
Markovian binary trees, which is supercritical for t G [0,1) and critical for 
t = 1, and let us denote with an additional subscript t the quantities defined 
above for this family of problems. Let us suppose that Rt := bt{e, •) + 6t(-,e) 
is irreducible for every t € [0, 1], and let p{JFy»j) be the spectral radius of the 



Jacobian of the Perron iteration as defined in (16). Then 
\im p{JFy^,t) -~ 



vfbi{yi,yi) w'^yi 



w'^hiljji,yi) vfyi 

where yi and are left and right Perron vectors of Ri. As a special case, if 
the vector w is a scalar multiple ofvi, the limit is 0. 

Proof. Let us define yt as the Perron vector of -ffj,*.* normalized so that \\yt\\i = 
1, and similarly v'[ as the left Perron vector of the same matrix, normalized so 
that \\yt\\i = 1. Since Hy»^t is irreducible, its left and right Perron vectors are 
analytical functions of t, and thus yt and Vt converge to yi and vi respectively. 
We have Hy*^i = Hq i = Ri. Notice that — > w'^{I — Ri), and that 

A\ = {b(e -y*,-) + b{-,e) - /)t ^ {R, - /)t. 
Moreover, since yi and vi span the right and left kernel of / — i?i , we have 

(I-i?i)(/-i?i)t=/-^. 

vi yi 

Additionally, we shall make use of the relation p{AB) — p{BA), valid for any A 
and B such that AB and BA are square matrices, in the first and second-to-last 
step of the following computation. 
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Putting all together, we get 



Therefore, we obtain 

\imp{JFy*^t) =P 



=P 



I- 



h{-,v*t ) - 




w'^btiyt,yt) 



bi{yi,yi)w'^ 
^w'^bi{yi,yi) 

' bi{yi,yi)w'^ 
^w'^bi{yi,yi 
1 



[Ri-iy I 



yivj 
vi yi 



1 - 



vibi{yi,yi) w^yi 
w'^bi{yi,yi) vfyi 




bi(jji,yi) 



□ 



For the normalization condition, the above result suggests taking w as the 
left Perron vector of 6(e, •) + 6(-,e). Indeed, this choice guarantees the local 
convergence of the Perron iteration for close-to-critical problems. Moreover we 
point out that, even though the convergence is linear, the speed of convergence 
increases as the MBT gets closer to critical; in particular, the convergence is 
superlinear in the critical case. 



6 Numerical experiments 

We compared the Perron iteration (PI) with the Newton method (NM) |3] and 
with the thicknesses algorithm (TH) [3]. As stated before, TH and PI are 
linearly convergent algorithms, while NM is a quadratically convergent one. All 
the experiments were performed using Matlab 7 (R14) on an Intel Xeon 2.80Ghz 
bi-processor. 

We applied the algorithms to the two test cases reported in [3]. The first 
one (El) is an MBT of size n = 9 depending on a parameter A, which is critical 
for A « 0.85 and supercritical for larger values of A. The second one (E2) is a 
MBT of size n — 3 depending on a parameter A, which is critical for A w 0.34 
and A « 0.84, and supercritical for the values inbetween. 

The only noteworthy issue in the implementation of PI is the method used 
for the computation of the maximal eigenvector. The classical methods are 



10 



Figure 1: CPU time in sec. (and number of iterations in brackets) for TH, NM 
and PI on several choices of A for El (top) and E2 (bottom) 



A 


TH 


NM 


PI 


0.86 


2.3935e+00 (11879) 


5.0932e-03 (14) 


4.9267e-03 (7) 


0.9 


6.5353e-01 (3005) 


4.2859e-03 (12) 


5.5756e-03 (8) 


1 


2.8049e-01 (1149) 


3.9009e-03 (11) 


5.5090e-03 (8) 


2 


9.2644e-02 (191) 


2.8453e-03 (8) 


5.5125e-03 (8) 


A 


TH 


NM 


PI 


0.5 


7.7003e-02 (132) 


2.3305e-03 (8) 


5.6983e-03 (11) 


0.7 


7.6503e-02 (135) 


2.1842e-03 (8) 


5.6081e-03 (11) 


0.8 


9.3603e-02 (313) 


2.4543e-03 (9) 


4.6166e-03 (9) 


0.84 


7.3060e-01 (4561) 


4.0001e-03 (13) 


4.1090e-03 (8) 



usually optimized for matrices of much larger size; however, here we deal with 
matrices of size n — 3 and n ~ 9, for which the complexity constants mat- 
ter. We compared several candidates (eigs, eig, the power method, a power 
method accelerated by repeated squaring of the matrix), and found that in our 
examples the fastest method to find the maximal eigenvector is computing the 
full eigenvector basis with [V, Lambda] =eig(P) and then selecting the maximal 
eigenvector. The picture should change for problems of larger size: eig takes 
O(n^) operations, while for instance eigs should take only O(n^) in typical 
cases. On the other hand, we point out that in absence of any structure (such 
as sparsity) in &(-,-)i forming the matrix &(w, •) or 6(-,w) for a new vector v, 
an operation which is required at every step in all known iterative algorithms, 
requires 0{n^) operations. Therefore, the CPU times are somehow indicative of 
the real complexity of the algorithms, but should be taken with a grain of salt. 

The stopping criterion was chosen to be \\x ~ a + b{x, x)\\ < ne, with e = 
10~^^, for all algorithms. 

The table in Figure[T]shows the results for several choices of A. The algorithm 
TH is clearly the slowest, taking far more CPU time than the two competitors. 
The different behavior of PI when approaching the critical cases is apparent: 
while the iterations for TH and NM increase, PI seems to be unaffected by 
the near-singularity of the problem, and in fact the iteration count decreases 
slightly. 

To show further results on the comparison between NM and PI, we report 
a number of graphs comparing the iteration count and CPU times of the two 
algorithms. The graphs are not cut at the critical values, but they extend to 
subcritical cases as well. It is an interesting point to note that when the MET 
is subcritical, and thus the minimal solution x* (extinction probability) is e, the 
two algorithms have a different behavior: NM (and TH as well) converges to e, 
while PI skips this solution and converges to a different solution x > e. This is 
because in the derivation of the Perron iteration we chose the solution a 7^ for 
the normalization equation, thus explicitly excluding the solution y = (i.e., 
X = e). 
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Figure [2] shows a plot of the iteration count of the two methods vs. different 
values of the parameter A. While in close-to-critical cases the iteration count for 
NM has a spike, the one for PI seems to decrease. However, the iteration count 
comparison is not fair since the steps of the two iterations require a different 
machine time. Figure |3] shows a similar plot, considering the CPU time instead 
of the iteration count. In order to achieve better accuracy, the plotted times are 
averages over 100 consecutive runs. 

The results now favor the Newton method in most experiments, but in close- 
to-critical cases the new method achieves better performance. The results are 
very close to each other, though, so it is to be expected that for larger input sizes 
or different implementations the differences in the performance of the eigensolver 
could lead to significant changes in the results. 

In order to highlight the performance difference in close-to-critical cases, we 
report in Figure |4] a plot with the CPU times sampled at a larger number of 
points around the most "interesting" regions of the previous graphs. 

The Jacobian ( [T6| had spectral radius less than 1 in all the above experi- 
ments, a condition which is needed to ensure the convergence of PI. However, 
this is not true for all possible MBTs. In fact, by setting the parameter A for 
El to much larger values, we encountered problematic cases in which PI did not 
converge. Specifically, starting from A ~ 78 the Jacobian ( 16 1 is larger than 1 
and PI does not converge. However, such cases are of little practical interest 
since they are highly supercritical MBTs, distant from the critical case, and 
thus they are easily solved with the traditional methods (NM or the customary 
functional iterations |2]) with a small number of iterations. 

The problem E2 is well-posed only for < A < 1, otherwise negative entries 
appear in 6, thus the above discussion does not apply. 

Along all the experiments reported above, all the matrices Hy appearing in 
the PI steps always turned out to have positive entries, even in the subcritical 
problems; thus their Perron vector and values were always well-defined and real. 



7 Conclusions and open issues 

We have proposed a new algorithm for solving the quadratic vector equation 
([T]), based on a Perron iteration. The algorithm performs well, both in terms 
of speed of convergence and accuracy, for close-to-critical problems where the 
classical methods are slower. 

Along the framework that we have exposed, several different choices are 
possible in the practical implementation of the new algorithm. 

One of them is the choice of the bilinear form •). Equation ^ and its 
solution depend only on the quadratic form b{t,t); however, there are different 
ways to extend it to a bilinear form b{s,t). This choice ultimately reflects a 
modeling aspect of the problem: when an individual spawns, it is transformed 
into two individuals in different states, and we may choose arbitrarily which of 
them is called the mother and which the child. 

As an example of how this choice affects the solution algorithms, changing 
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Figure 4: Detailed views from Figure |3j CPU time ( in sec.) vs. parameter A 
for El (top) and E2 (middle and bottom) 
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the bilinear form may transform the depth algorithm into the order one and vice 
versa. The algorithms we proposed depend on the actual choice of the bilinear 
extension of the quadratic form b{t, t), and the convergence speed is affected by 
this decision. 

A second choice is the normalization of the computed Perron vector: different 
approaches may be attempted — for instance, minimization of the 1-norm, of 
the 2-norm, or orthogonality of the residual of ^ with respect to a suitably 
chosen vector — although it is not clear whether we can improve the results of 
the normalization presented here. 

A third choice, crucial in the computational experiments, is the method used 
to compute the Perron vector. For moderate sizes of the problem, it is cheaper 
to do a full eigendecomposition of the matrix and extract the eigenvalue with 
maximum modulus, but for larger problems it pays off to use different specific 
methods for its computation. 

All these variants deserve to be better understood, and are now under our 
investigation. 
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